%fan-beam geometry parameter setting %
nx=128;ny=nx;% number of pixels: nx*ny
dx=single(250/nx);dy=dx;% pixel size: dx*dy
nv=15;% number of views (projections)
SO=single(1000);% the distance from source to isocenter
OD=single(500);% the distance from isocenter to detector
nd=512;% number of detector
dy_det=single(0.388*1024/nd);% size of detector
fan_angle = atan(512/1500);%fan angle
sd_phi=single(1*(pi+fan_angle)/nv*(0:nv-1));% view angles
y_os=single(0*dy_det);% isocenter offset with respect to the detector center
y_det=single(((-nd/2:nd/2-1)+0.5)*dy_det)+y_os;% detector coordinate
y_det2=single((-nd/2:nd/2)*dy_det)+y_os;
% normalize so that dx=dy=1 %
scale=dx;
SO=SO/scale;OD=OD/scale;y_det=y_det/scale;
dy_det=dy_det/scale;y_os=y_os/scale;





%set Parameter
nt=1;
para =  setPara(SO,OD,scale,nx,ny,nv,nt,dy_det,y_os,sd_phi,y_det,y_det2);

%

% object%
x0 = phantom(nx);
x0(x0<0)=0;
x0 = single(x0);
X0 = x0(:);

% X-ray Transform %
disp( 'X-ray Transform');
para.version=uint32(0);para.GPU=uint32(0); % Siddon's algorithm - CPU
tic;y_siddon_cpu=Ax_fan_mf(X0,para);toc;

%iterative recon;
x0 = single(zeros(size(X0)));
prjimg = y_siddon_cpu;